clear
rep=1;

%% Problem parameters
d=2;
x_bound=repmat([-6,6],[d,1]);
precision = 0.5;
NewRegionNum = 2;
TotalBudget = 1000;
MaxPartitionDepth_d = ceil(log((x_bound(:,2)-x_bound(:,1))./precision)...
    ./log(NewRegionNum));
MaxPartitionDepth = sum(MaxPartitionDepth_d);

Problem.Dimension = d;
Problem.Domain = reshape(x_bound',[1,d*2]);
Problem.Partition = @(Region,SampleSet)Partition_NP_BoxConstraint(Region, NewRegionNum, MaxPartitionDepth, SampleSet);
Problem.Sampling = @Sampling_BoxConstraint;

%% Algorithm parameters
AlgorithmM.n0 = 4;
AlgorithmM.QuantileLevel = 0.18;
AlgorithmM.NewBudget = 3;
AlgorithmM.MaxSampleSize = 15;
AlgorithmM.StopCriteria = [1,TotalBudget];

ClearRadial0 = (x_bound(:,2)-x_bound(:,1))./(NewRegionNum.^(MaxPartitionDepth_d));
ClearRadial = sqrt(sum(ClearRadial0.^2)); % the maximum distance within a non-partitionable region

%% STAGE ONE
time1=tic;
for i=1:rep
[RegionSet, SampleSet, RegionSampleId] = PRS_MMO( Problem, AlgorithmM );
[Optima] = Niching(RegionSet,SampleSet,RegionSampleId,ClearRadial);
end
toc(time1)

%% FIGURES
if rep == 1
    if d == 2
        figure()
        hold on
        rectangle('Position',[Problem.Domain([1,3]),Problem.Domain([2,4])-Problem.Domain([1,3])])
        scatter(SampleSet(:,2),SampleSet(:,3),2,SampleSet(:,1),'filled');
        for i=1:length(RegionSet)
            rectangle('Position',[RegionSet(i,[3,5]),RegionSet(i,[4,6])-RegionSet(i,[3,5])])
        end
        title('','Interpreter','latex','String',['$N=',num2str(TotalBudget),'$'])
        xlabel('','Interpreter','latex','String','$x_1$')
        ylabel('','Interpreter','latex','String','$x_2$')
        set(gca,'Fontname','Times New Roman','FontSize',16);
        for i=1:size(Optima,1)
            plot(Optima(i,2),Optima(i,3),'x');
        end
        hold off
    end
end

%% STAGE 2 LOCAL SEARCH
Optima_Stage2 = Optima;
TotalBudget_Stage2 = zeros(size(Optima,1),1);
Path = cell(size(Optima,1),1);
for i = 1:size(Optima,1)
    step = 0.005;
    currentbest = Optima(i,:);
    Path{i} = currentbest;
    while 1
        neighbor = zeros(2*d,d);
        neighbor(1:(2*d+2):(2*d*d)) = 1;
        neighbor(2:(2*d+2):(2*d*d)) = -1;
%         neighbor = (fullfact(ones(1,d).*2)-1.5).*2;
        NumNeighbor = size(neighbor,1);
        neighbor = neighbor.*step;
        neighbor = neighbor + repmat(currentbest(2:end),[NumNeighbor,1]);
        fvalue = ObjValue(neighbor);
        TotalBudget_Stage2(i) = TotalBudget_Stage2(i)+NumNeighbor;
        [~,best_temp] = min([fvalue;currentbest(1)]);
        if best_temp==(NumNeighbor+1)
            break
        else
            currentbest = [fvalue(best_temp),neighbor(best_temp,:)];
            Path{i} = [Path{i};currentbest];
        end
    end
    Optima_Stage2(i,:) = currentbest;
end

%% FIGURES
figure()
hold on
rectangle('Position',[Problem.Domain([1,3]),Problem.Domain([2,4])-Problem.Domain([1,3])])
for i=1:size(Optima_Stage2,1)
    plot(Path{i}(:,2), Path{i}(:,3), '-.')
    plot(Path{i}(end,2), Path{i}(end,3), 'o')
end
hold off
xlabel('','Interpreter','latex','String','$x_1$')
ylabel('','Interpreter','latex','String','$x_2$')
title('Local search')
set(gca,'Fontname','Times New Roman','FontSize',16);
%%
% save('data.mat')

%% RANDOMLY SAMPLING
SampleSet_random = Problem.Sampling(TotalBudget,Problem.Domain);
RegionSampleId_random = cell(1,1);
RegionSampleId_random{1} = (1:TotalBudget)';

%% 660 new budget are needed
ClearRadial = 0.86;%0.86-3.8
[Optima_random] = Niching([1,0,Problem.Domain],SampleSet_random,RegionSampleId_random,ClearRadial);

figure()
hold on
rectangle('Position',[Problem.Domain([1,3]),Problem.Domain([2,4])-Problem.Domain([1,3])])
scatter(SampleSet_random(:,2),SampleSet_random(:,3),2,SampleSet_random(:,1),'filled');
for i=1:size(Optima_random,1)
    plot(Optima_random(i,2),Optima_random(i,3),'x');
end
hold off